In this script we will provide some intuition for harmonic regression, and walk through how to reconstruct temporally dynamic coefficients when fitting step selection (or other) models with harmonic terms. We will use the example of a model fitted with two pairs of harmonics. We will then use the coefficients to reconstruct the selection surfaces when we interact linear and quadratic terms with the harmonics.

Load packages


options(scipen=999)

library(tidyverse)
packages <- c("lubridate", "tictoc", "TwoStepCLogit", "beepr", "scales")
walk(packages, require, character.only = T)

Harmonics are sine and cosine functions are periodic functions that repeat once per cycle, \(T\). When combined, they can be used to model complex periodic patterns, and with many harmonics can create very flexible functions.

Harmonics can be considered similarly to splines or basis functions in that they are additive - if we add them all together we get a single, flexible function.

Let’s say we have a single pair of harmonics, \(s_1 = \sin(2\pi t / T)\), and \(c_1 = \cos(2\pi t / T)\), where \(t \in T\). For a yearly cycle for instance, \(t\) would be indexed as the day of the year, and \(T = 365\), although \(t\) does not need to be integer-valued, and can be arbitrarily fine. In the case of the model we are fitting here, \(t\) is the hour of the day, for \(T\) is the number of hours in a day, 24.


# create a sequence of hours
t <- seq(0,23,1)
T <- 24

# create the harmonic terms
sin1 <- sin(2 * pi * t / T)
cos1 <- cos(2 * pi * t / T)

# plot the harmonic terms
plot(t, sin1, type = "l", col = "red", 
     xlab = "Hour of the day", ylab = "Value", 
     main = "Harmonic terms", ylim = c(-2.5,2.5))
lines(t, cos1, col = "blue")
lines(t, rep(0,T), lty = "dashed")

Now what we can do is add them together to create a single function.


# create a single function
f_harmonic <- sin1 + cos1

# plot the function
plot(t, f_harmonic, type = "l", col = "black", 
     xlab = "Hour of the day", ylab = "Value", 
     main = "Harmonic function", 
     ylim = c(-2.5,2.5))
lines(t, rep(0,T), lty = "dashed")

It looks similar, although now the amplitude is greater (up to about 1.4 now), and that the peak has shifted to be an average of the sine and cosine functions (peak is now at 3 ( = 0 + 6 / 2). We can shift the location of the peak by giving one of the terms more influence.


# create a single function
f_harmonic <- 2*sin1 + cos1

# plot the function
plot(t, f_harmonic, type = "l", col = "black", 
     xlab = "Hour of the day", ylab = "Value", 
     main = "Harmonic function", ylim = c(-2.5,2.5))
lines(t, rep(0,T), lty = "dashed")

If we give more weight to the sine function, the peak will shift to the right, which is closer to the peak of the sine function, and vice versa. We can see it has also changed the amplitude, which now peaks above y = 2.

We can also add a constant term, which will just shift the entire function up or down.


# create a single function
f_harmonic <- 1 + 2*sin1 + cos1

# plot the function
plot(t, f_harmonic, type = "l", col = "black", 
     xlab = "Hour of the day", ylab = "Value", 
     main = "Harmonic function", ylim = c(-2.5,2.5))
lines(t, rep(0,T), lty = "dashed")

We can see that adding a constant, and by weighting the sin and cosine terms we can start to create a flexible function. If we only use 1 pair of harmonics we can only have one period, but if we use more pairs we can have more periods, and therefore more flexibility. Let’s add some more pairs of harmonics.


# create two more harmonic terms
sin2 <- sin(4 * pi * t / T)
cos2 <- cos(4 * pi * t / T)

# plot the harmonic terms
plot(t, sin2, type = "l", col = "red", 
     xlab = "Hour of the day", ylab = "Value", 
     main = "Harmonic terms", ylim = c(-2.5,2.5))
lines(t, cos2, col = "blue")
lines(t, rep(0,T), lty = "dashed")

Again we can add these together, and weight to shift the peak and change the amplitude.


# create a single function
f_harmonic <- 0.5*sin2 + 1.5*cos2

# plot the function
plot(t, f_harmonic, type = "l", col = "black", 
     xlab = "Hour of the day", ylab = "Value", 
     main = "Harmonic function", ylim = c(-2.5,2.5))
lines(t, rep(0,T), lty = "dashed")

Now we can start to add all of our components together to create flexible functions (and there is really no limit to this - infinitely flexible functions can be created with infinitely many harmonic terms, which is called a Fourier series - for some very cool intuition about Fourier series check out: https://www.youtube.com/watch?v=r6sGWTCMz2k).


# create a single function
f_harmonic <- 1 + 0.5*sin1 + 2*cos1 + 0.75*sin2 + 1.25*cos2

# plot the function
plot(t, f_harmonic, type = "l", col = "black", 
     xlab = "Hour of the day", ylab = "Value", 
     main = "Harmonic function")
lines(t, rep(0,T), lty = "dashed")

If you play around with the weights you can estimate some really funky functions, and the function will just get more flexible the more harmonic terms that you add.


# create more harmonic terms...
sin3 <- sin(6 * pi * t / T)
cos3 <- cos(6 * pi * t / T)
sin4 <- sin(8 * pi * t / T)
cos4 <- cos(8 * pi * t / T)

# create a single function
f_harmonic <- 0.5 + 0.5*sin1 + 1.5*cos1 + 2*sin2 + 1.25*cos2 + 0.25*sin3 + 0.5*cos3 + 0.1*sin4 + 2*cos4

# plot the function
plot(t, f_harmonic, type = "l", col = "black", 
     xlab = "Hour of the day", ylab = "Value", 
     main = "Harmonic function")
lines(t, rep(0,T), lty = "dashed")

Hopefully now you have some intuition about what we are doing with the harmonic terms in our step selection function, in that we are estimating the weights of the sine and cosine functions, which are the coefficients from our fitted models. Pretty cool!

Let’s have a look at some real data now, and fit a model with the harmonics.

Buffalo example

Importing buffalo data for model fitting, and creating harmonic terms using the time at the end of the step.


buffalo_data_all <- read_csv(
  "outputs/buffalo_popn_GvM_covs_ST_KDEmem1000_allOPTIM_10rs_2024-02-05.csv")
## Rows: 1082697 Columns: 26
## ── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: ","
## dbl  (22): id, burst_, x1_, x2_, y1_, y2_, sl_, ta_, dt_, hour_t2, step_id_, y, ndvi_temporal, veg_woody, veg_herby...
## lgl   (1): case_
## dttm  (3): t1_, t2_, t2_rounded
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

buffalo_data_all <- buffalo_data_all %>% 
  mutate(t1_ = lubridate::with_tz(buffalo_data_all$t1_, tzone = "Australia/Darwin"),
         t2_ = lubridate::with_tz(buffalo_data_all$t2_, tzone = "Australia/Darwin"))

buffalo_data_all <- buffalo_data_all %>%
  mutate(id_num = as.numeric(factor(id)),
         step_id = step_id_,
         x1 = x1_, x2 = x2_,
         y1 = y1_, y2 = y2_,
         t1 = t1_,
         t1_rounded = round_date(buffalo_data_all$t1_, "hour"),
         hour_t1 = hour(t1_rounded),
         t2 = t2_,
         t2_rounded = round_date(buffalo_data_all$t2_, "hour"),
         hour_t2 = hour(t2_rounded),
         hour = hour_t2,
         yday = yday(t1_),
         year = year(t1_),
         month = month(t1_),
         sl = sl_,
         log_sl = log(sl_),
         ta = ta_,
         cos_ta = cos(ta_),
         # scale canopy cover from 0 to 1
         canopy_01 = canopy_cover/100,
         # recover the natural scale of previous space use density
         kde_memory_density = exp(kde_memory_density_log),
         
         # here we create the harmonic terms for the hour of the day
         # for seasonal effects, change hour to day of the year (which is t), 
         # and 24 to 365 (which is T)
         hour_s1 = sin(2*pi*hour/24),
         hour_s2 = sin(4*pi*hour/24),
         hour_s3 = sin(6*pi*hour/24),
         hour_c1 = cos(2*pi*hour/24),
         hour_c2 = cos(4*pi*hour/24),
         hour_c3 = cos(6*pi*hour/24))

# to just select a single year of data
buffalo_data_all <- buffalo_data_all %>% filter(t1 < "2019-07-25 09:32:42 ACST")
buffalo_ids <- unique(buffalo_data_all$id)

Fitting the model with two pairs of harmonics

Creating a data matrix

First we create a data matrix to be provided to the model, and then we scale and centre the full data matrix, with respect to each of the columns. That means that all variables are scaled and centred after the data has been split into wet and dry season data, and also after creating the quadratic and harmonic terms (when using them).

This is where we interact the harmonic terms with our covariates. As we have already created the harmonic terms for the hour of the day (s1, c1, s2, etc), we just interact (multiply) these with each of the covariates, including the quadratic terms and movement parameters. As before, we create the data matrix with all quadratic and harmonic terms, and then scale the matrix by each column, and store the scaling and centering variables to reconstruct the natural scale coefficients (i.e. if we were to fit the model without centering and scaling).


months_wet <- c(1:4, 11:12)
buffalo_ids <- unique(buffalo_data_all$id)

# buffalo_data <- buffalo_data_all %>% filter(month %in% months_wet) # wet season
buffalo_data <- buffalo_data_all %>% filter(!month %in% months_wet) # dry season

buffalo_data_matrix_unscaled <- buffalo_data %>% transmute(
  
  ndvi = ndvi_temporal,
  ndvi_s1 = ndvi_temporal * hour_s1,
  ndvi_s2 = ndvi_temporal * hour_s2,
  ndvi_c1 = ndvi_temporal * hour_c1,
  ndvi_c2 = ndvi_temporal * hour_c2,
  
  ndvi_sq = ndvi_temporal ^ 2,
  ndvi_sq_s1 = (ndvi_temporal ^ 2) * hour_s1,
  ndvi_sq_s2 = (ndvi_temporal ^ 2) * hour_s2,
  ndvi_sq_c1 = (ndvi_temporal ^ 2) * hour_c1,
  ndvi_sq_c2 = (ndvi_temporal ^ 2) * hour_c2,
  
  canopy = canopy_01,
  canopy_s1 = canopy_01 * hour_s1,
  canopy_s2 = canopy_01 * hour_s2,
  canopy_c1 = canopy_01 * hour_c1,
  canopy_c2 = canopy_01 * hour_c2,
  
  canopy_sq = canopy_01 ^ 2,
  canopy_sq_s1 = (canopy_01 ^ 2) * hour_s1,
  canopy_sq_s2 = (canopy_01 ^ 2) * hour_s2,
  canopy_sq_c1 = (canopy_01 ^ 2) * hour_c1,
  canopy_sq_c2 = (canopy_01 ^ 2) * hour_c2,
  
  slope = slope,
  slope_s1 = slope * hour_s1,
  slope_s2 = slope * hour_s2,
  slope_c1 = slope * hour_c1,
  slope_c2 = slope * hour_c2,
  
  herby = veg_herby,
  herby_s1 = veg_herby * hour_s1,
  herby_s2 = veg_herby * hour_s2,
  herby_c1 = veg_herby * hour_c1,
  herby_c2 = veg_herby * hour_c2,
  
  spatial_memory = kde_memory_density,
  spatial_memory_s1 = kde_memory_density * hour_s1,
  spatial_memory_s2 = kde_memory_density * hour_s2,
  spatial_memory_c1 = kde_memory_density * hour_c1,
  spatial_memory_c2 = kde_memory_density * hour_c2,

  spatial_memory_sq = kde_memory_density ^ 2,
  spatial_memory_sq_s1 = (kde_memory_density ^ 2) * hour_s1,
  spatial_memory_sq_s2 = (kde_memory_density ^ 2) * hour_s2,
  spatial_memory_sq_c1 = (kde_memory_density ^ 2) * hour_c1,
  spatial_memory_sq_c2 = (kde_memory_density ^ 2) * hour_c2,
  
  step_l = sl,
  step_l_s1 = sl * hour_s1,
  step_l_s2 = sl * hour_s2,
  step_l_c1 = sl * hour_c1,
  step_l_c2 = sl * hour_c2,

  log_step_l = log_sl,
  log_step_l_s1 = log_sl * hour_s1,
  log_step_l_s2 = log_sl * hour_s2,
  log_step_l_c1 = log_sl * hour_c1,
  log_step_l_c2 = log_sl * hour_c2,

  cos_turn_a = cos_ta,
  cos_turn_a_s1 = cos_ta * hour_s1,
  cos_turn_a_s2 = cos_ta * hour_s2,
  cos_turn_a_c1 = cos_ta * hour_c1,
  cos_turn_a_c2 = cos_ta * hour_c2)

buffalo_data_matrix_scaled <- scale(buffalo_data_matrix_unscaled)

# store the scaling and centering attribues to recover the natural scale coefficients
mean_vals <- attr(buffalo_data_matrix_scaled, "scaled:center")
sd_vals <- attr(buffalo_data_matrix_scaled, "scaled:scale")
scaling_attributes <- data.frame(variable = names(buffalo_data_matrix_unscaled), 
                                 mean = mean_vals, sd = sd_vals)

buffalo_data_scaled_2p <- data.frame(id = buffalo_data$id,  
                                     step_id = buffalo_data$step_id, 
                                     y = buffalo_data$y, 
                                     buffalo_data_matrix_scaled)

Formula with two pairs of harmonics

We add in all the terms in the model formula. We can see that we are adding all of these together, which allows us to add them together later to reconstruct the harmonic function.


formula_twostep <- y ~ 
  
  ndvi +
  ndvi_s1 +
  ndvi_s2 +
  ndvi_c1 +
  ndvi_c2 +
  
  ndvi_sq +
  ndvi_sq_s1 +
  ndvi_sq_s2 +
  ndvi_sq_c1 +
  ndvi_sq_c2 +
  
  canopy +
  canopy_s1 +
  canopy_s2 +
  canopy_c1 +
  canopy_c2 +
  
  canopy_sq +
  canopy_sq_s1 +
  canopy_sq_s2 +
  canopy_sq_c1 +
  canopy_sq_c2 +
  
  slope +
  slope_s1 +
  slope_s2 +
  slope_c1 +
  slope_c2 +
  
  herby +
  herby_s1 +
  herby_s2 +
  herby_c1 +
  herby_c2 +

  spatial_memory +
  spatial_memory_s1 +
  spatial_memory_s2 +
  spatial_memory_c1 +
  spatial_memory_c2 +

  spatial_memory_sq +
  spatial_memory_sq_s1 +
  spatial_memory_sq_s2 +
  spatial_memory_sq_c1 +
  spatial_memory_sq_c2 +
  
  step_l +
  step_l_s1 +
  step_l_s2 +
  step_l_c1 +
  step_l_c2 +
  
  log_step_l +
  log_step_l_s1 +
  log_step_l_s2 +
  log_step_l_c1 +
  log_step_l_c2 +

  cos_turn_a +
  cos_turn_a_s1 +
  cos_turn_a_s2 +
  cos_turn_a_c1 +
  cos_turn_a_c2 +
  
  strata(step_id) +
  cluster(id)

Fitting the model

As we have already fitted the model, we will load it here, but if the model_fit object file doesn’t exist, it will run the model fitting code. Be careful here that if you change the model formula, you will need to delete or rename the model_fit file to re-run the model fitting code, otherwise it will just load the previous model.


if(file.exists("outputs/model_twostep_2p_harms_dry.rds")) {
  
  model_twostep_2p_harms <- readRDS("outputs/model_twostep_2p_harms_dry.rds")
  
} else {

  tic()
  model_twostep_2p_harms <- Ts.estim(formula = formula_twostep,
           data = buffalo_data_scaled_2p,
           all.m.1 = TRUE,
           D = "UN(1)",
           itermax = 10000)
  toc()
  
  # save model object
  saveRDS(model_twostep_2p_harms, file = "outputs/model_twostep_2p_harms_dry.rds")
  
  beep(sound = 2)

}

Check the model output

It will output a massive covariance matrix, which can be ignored for our purposes at the moment. We are interested in the coefficients, which are stored in the beta attribute of the model object.


model_twostep_2p_harms$beta
##                 ndvi              ndvi_s1              ndvi_s2              ndvi_c1              ndvi_c2 
##           1.18745765           0.39442696           0.01121123          -0.95521148           1.25104461 
##              ndvi_sq           ndvi_sq_s1           ndvi_sq_s2           ndvi_sq_c1           ndvi_sq_c2 
##          -1.33022674          -0.24830951          -0.01598593           0.12688921          -0.46560189 
##               canopy            canopy_s1            canopy_s2            canopy_c1            canopy_c2 
##           0.27443845           0.16425934          -0.14215576          -0.03190371           0.25099088 
##            canopy_sq         canopy_sq_s1         canopy_sq_s2         canopy_sq_c1         canopy_sq_c2 
##          -0.34884092          -0.09769922           0.04173728          -0.16166660          -0.16650068 
##                slope             slope_s1             slope_s2             slope_c1             slope_c2 
##          -0.16917086          -0.01981839          -0.05097721          -0.03595167          -0.05467170 
##                herby             herby_s1             herby_s2             herby_c1             herby_c2 
##          -0.02374414          -0.02758919           0.01137859           0.13312922          -0.04304990 
##       spatial_memory    spatial_memory_s1    spatial_memory_s2    spatial_memory_c1    spatial_memory_c2 
##           2.29925142          -0.48174644          -0.69537479          -0.61650328           0.48955057 
##    spatial_memory_sq spatial_memory_sq_s1 spatial_memory_sq_s2 spatial_memory_sq_c1 spatial_memory_sq_c2 
##          -1.39487590           0.17180133           0.55104843           0.60512613          -0.46486791 
##               step_l            step_l_s1            step_l_s2            step_l_c1            step_l_c2 
##          -0.34761318          -0.10567643          -0.61073611          -0.15348490          -0.32241309 
##           log_step_l        log_step_l_s1        log_step_l_s2        log_step_l_c1        log_step_l_c2 
##           0.23121452          -0.13842027          -0.13823993          -0.23188181           0.06039618 
##           cos_turn_a        cos_turn_a_s1        cos_turn_a_s2        cos_turn_a_c1        cos_turn_a_c2 
##          -0.01128474          -0.03790357          -0.18570605          -0.07010761          -0.02699471

Create a dataframe of the coefficients with the scaling attributes, and return the coefficients to their natural scale by dividing by the scaling factor (standard deviation).

As we can see, we have a coefficient for each covariate by itself, and then one for each of the harmonics. These are the ‘weights’ that we played around with above, and we reconstruct them in exactly the same way. We also have the coefficients for the quadratic terms and the interactions with the harmonics, which I have denoted as ndvi_sq for instance. We will come back to these when we look at the selection surfaces.


# creating data frame of model coefficients
coefs_clr <- data.frame(coefs = names(model_twostep_2p_harms$beta), 
                        value = model_twostep_2p_harms$beta)

coefs_clr$scale_sd <- scaling_attributes$sd
coefs_clr <- coefs_clr %>% mutate(value_nat = value / scale_sd)
head(coefs_clr, 10)
##                 coefs       value  scale_sd    value_nat
## ndvi             ndvi  1.18745765 0.1429054   8.30939941
## ndvi_s1       ndvi_s1  0.39442696 0.2371431   1.66324473
## ndvi_s2       ndvi_s2  0.01121123 0.2372765   0.04724962
## ndvi_c1       ndvi_c1 -0.95521148 0.2426054  -3.93730465
## ndvi_c2       ndvi_c2  1.25104461 0.2425239   5.15843931
## ndvi_sq       ndvi_sq -1.33022674 0.1019530 -13.04744654
## ndvi_sq_s1 ndvi_sq_s1 -0.24830951 0.1065570  -2.33029697
## ndvi_sq_s2 ndvi_sq_s2 -0.01598593 0.1076953  -0.14843668
## ndvi_sq_c1 ndvi_sq_c1  0.12688921 0.1105880   1.14740494
## ndvi_sq_c2 ndvi_sq_c2 -0.46560189 0.1097234  -4.24341508

Let’s have a look at herbaceous vegetation as it didn’t have a quadratic term.

Firstly, as above, we need a sequence of values that covers a full period (or the period that we want to construct the function over, which can be more or less than 1 period). The sequence can be arbitrarily finely spaced. The smaller the increment the smoother the function will be for plotting. When simulating data from the temporally dynamic coefficients, however, one would use the increment that relates to the data collection and model fitting (i.e. one hour in this case).

Here we’ll use a finer resolution than above for smoother plotting


hour_seq <- seq(0,23.9,0.1)

Now we can reconstruct the harmonic function using the formula that we put into our model by interacting the harmonic terms with each of the covariates, for a single covariate, let’s say herbaceous vegetation, this would be written down as:

\[ f = \beta_{herby} + \beta_{herby\_s1} \sin\left(\frac{2\pi t}{24}\right) + \beta_{herby\_c1} \cos\left(\frac{2\pi t}{24}\right) + \beta_{herby\_s2} \sin\left(\frac{4\pi t}{24}\right) + \beta_{herby\_c2} \cos\left(\frac{4\pi t}{24}\right) \]

We can precompute the harmonic terms to make it a bit neater (as we did for the design matrix before fitting our model)


t <- hour_seq
T <- 24

sin1 <- sin(2 * pi * t / T)
cos1 <- cos(2 * pi * t / T)
sin2 <- sin(4 * pi * t / T)
cos2 <- cos(4 * pi * t / T)

Now we just pull out the relevant coefficients (which are scalars (single numbers), like the weights we used above), multiply them by the harmonic terms (which are vectors representing the function through time), and add them together.


herby_harmonic_function <- coefs_clr$value[which(coefs_clr$coefs == "herby")] +
  coefs_clr$value[which(coefs_clr$coefs == "herby_s1")] * sin1 +
  coefs_clr$value[which(coefs_clr$coefs == "herby_c1")] * cos1 +
  coefs_clr$value[which(coefs_clr$coefs == "herby_s2")] * sin2 +
  coefs_clr$value[which(coefs_clr$coefs == "herby_c2")] * cos2

# plot the function
plot(t, herby_harmonic_function, type = "l", col = "black", 
     xlab = "Hour of the day", ylab = "Coefficient value (scaled)", 
     main = "Herbaceous vegetation")
lines(t, rep(0,length(t)), lty = "dashed")

From these we can see that buffalo select for herbaceous vegetation in the early morning and late afternoon, but have a reasonably strong avoidance of it in the middle of the day (and therefore a selection for woody vegetation). This aligns with what we know about buffalo behaviour (and the climate in Northern Australia’s tropical savannas), in that they are likely to be seeking shelter during this time due to high temperatures and sun.

Reconstructing coefficients using matrix multiplication

Now that we know how the coefficients of the harmonics can be combined to form a temporally dynamic function, we can use a shortcut, which is matrix multiplication.

First we create a matrix of the harmonics, which is just the sin and cos terms for each harmonic, and then we can multiply this by the coefficients to get the function. When we use two pairs of harmonics we will have 5 coefficients for each covariate (linear + 2 sin and 2 cos), so there will be 5 columns in the matrix.

For matrix multiplication, the number of columns in the first matrix must be equal to the number of rows in the second matrix. The result will then have the same number of rows as the first matrix and the same number of columns as the second matrix.

Or in other words, if we have a 24 x 5 matrix of harmonics and a 5 x 1 matrix of coefficients, we will get a 24 x 1 matrix of the function, which corresponds to our 24 hours of the day (and the middle numbers be the same).


# we'll return the increments back to 1 hour for this example, 
# just so everything is a bit clearer
hour_seq <- seq(0,23,1)

hour_harmonics_matrix <- as.matrix(data.frame("linear_term" = rep(1, length(hour_seq)),
                                              "hour_s1" = sin(2*pi*hour_seq/24),
                                              "hour_s2" = sin(4*pi*hour_seq/24),
                                              "hour_c1" = cos(2*pi*hour_seq/24),
                                              "hour_c2" = cos(4*pi*hour_seq/24)))

# now have a 24 x 5 matrix
head(hour_harmonics_matrix)
##      linear_term   hour_s1   hour_s2   hour_c1                    hour_c2
## [1,]           1 0.0000000 0.0000000 1.0000000  1.00000000000000000000000
## [2,]           1 0.2588190 0.5000000 0.9659258  0.86602540378443870761060
## [3,]           1 0.5000000 0.8660254 0.8660254  0.50000000000000011102230
## [4,]           1 0.7071068 1.0000000 0.7071068  0.00000000000000006123032
## [5,]           1 0.8660254 0.8660254 0.5000000 -0.49999999999999977795540
## [6,]           1 0.9659258 0.5000000 0.2588190 -0.86602540378443870761060

We can now pull out the coefficients for each covariate and multiply by the matrix to get our temporally dynamic function. I’ll do this bit by bit initially.

Starting with the scaled NDVI covariate


# this subsets the coefficients for the covariate of interest 
# (grepl uses string matching, and we include ndvi but exclude the quadratic term)
ndvi_coefs <- coefs_clr %>% filter(grepl("ndvi", coefs) & !grepl("sq", coefs)) 
ndvi_coefs
##           coefs       value  scale_sd   value_nat
## ndvi       ndvi  1.18745765 0.1429054  8.30939941
## ndvi_s1 ndvi_s1  0.39442696 0.2371431  1.66324473
## ndvi_s2 ndvi_s2  0.01121123 0.2372765  0.04724962
## ndvi_c1 ndvi_c1 -0.95521148 0.2426054 -3.93730465
## ndvi_c2 ndvi_c2  1.25104461 0.2425239  5.15843931
       
# pull out the coefficients - this comes out as a column vector, 
# which we'll convert to a matrix
ndvi_coefs_matrix <- matrix(ndvi_coefs %>% pull(value))
ndvi_coefs_matrix
##             [,1]
## [1,]  1.18745765
## [2,]  0.39442696
## [3,]  0.01121123
## [4,] -0.95521148
## [5,]  1.25104461

# now we can multiply the matrix by the coefficients to get the function
ndvi_harmonic_function <- hour_harmonics_matrix %*% ndvi_coefs_matrix
# ndvi_harmonic_function

# plot the function
plot(hour_seq, ndvi_harmonic_function, type = "l", col = "black", 
     xlab = "Hour of the day", ylab = "Coefficient value", 
     main = "NDVI - linear term")
lines(hour_seq, rep(0,length(hour_seq)), lty = "dashed")

Now we can just repeat the process for all of the covariates, and put it all into a dataframe.


# back to the finer resolution for smoother plotting
hour <- seq(0,23.9,0.1)

hour_harmonics_matrix <- as.matrix(data.frame("linear_term" = rep(1, length(hour)),
                                              "hour_s1" = sin(2*pi*hour/24),
                                              "hour_s2" = sin(4*pi*hour/24),
                                              "hour_c1" = cos(2*pi*hour/24),
                                              "hour_c2" = cos(4*pi*hour/24)))

harmonics_df_2p <- data.frame(
  "hour" = hour,
  "ndvi" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("ndvi", coefs) & !grepl("sq", coefs)) %>% 
             pull(value)),
  "ndvi_2" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("ndvi_sq", coefs)) %>% pull(value)),
  "canopy" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("canopy", coefs) & !grepl("sq", coefs)) %>% 
             pull(value)),
  "canopy_2" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("canopy_sq", coefs)) %>% pull(value)),
  "slope" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("slope", coefs)) %>% pull(value)),
  "herby" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("herby", coefs)) %>% pull(value)),
  "memory" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("memory", coefs) & !grepl("sq", coefs)) %>% 
             pull(value)),
  "memory_2" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("memory_sq", coefs)) %>% pull(value)),
  "sl" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("step", coefs) & !grepl("log", coefs)) %>% 
             pull(value)),
  "log_sl" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("log_step", coefs)) %>% pull(value)),
  "cos_ta" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("turn_a", coefs)) %>% pull(value))
  )

head(harmonics_df_2p)
##   hour     ndvi    ndvi_2    canopy   canopy_2      slope      herby   memory  memory_2         sl     log_sl
## 1  0.0 1.483291 -1.668939 0.4935256 -0.6770082 -0.2597942 0.06633518 2.172299 -1.254618 -0.8235112 0.05972889
## 2  0.1 1.492815 -1.675681 0.4900525 -0.6770977 -0.2628937 0.06622186 2.122835 -1.220851 -0.8577465 0.04886724
## 3  0.2 1.499561 -1.681229 0.4859317 -0.6766250 -0.2658112 0.06613404 2.072563 -1.186309 -0.8909046 0.03802145
## 4  0.3 1.503533 -1.685581 0.4811834 -0.6755962 -0.2685398 0.06607018 2.021594 -1.151077 -0.9229000 0.02721430
## 5  0.4 1.504745 -1.688739 0.4758296 -0.6740188 -0.2710729 0.06602851 1.970045 -1.115247 -0.9536502 0.01646880
## 6  0.5 1.503219 -1.690708 0.4698938 -0.6719016 -0.2734045 0.06600701 1.918035 -1.078908 -0.9830758 0.00580812
##       cos_ta
## 1 -0.1083871
## 2 -0.1190373
## 3 -0.1295384
## 4 -0.1398633
## 5 -0.1499856
## 6 -0.1598791

# turning into a long data frame for plotting with ggplot
harmonics_df_2p_long <- pivot_longer(harmonics_df_2p, cols = !1, names_to = "coef")

Now we can plot all these together to see the temporal dynamics of the coefficients. Keep in mind that these are the scaled coefficients.


ggplot() +
  geom_path(data = harmonics_df_2p_long,
            aes(x = hour, y = value, colour = coef)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_y_continuous(expression(Time-varying~parameter~values~beta)) +
  scale_x_continuous("Hour", breaks = seq(0,24,2)) +
  scale_color_discrete("Estimate") +
  theme_classic() +
  theme(legend.position = "bottom")

We can do the same with the natural scale coefficients.


harmonics_nat_df_2p <- data.frame(
  "hour" = hour,
  "ndvi" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("ndvi", coefs) & !grepl("sq", coefs)) %>% 
             pull(value_nat)),
  "ndvi_2" = hour_harmonics_matrix %*%
    matrix(coefs_clr %>% filter(grepl("ndvi_sq", coefs)) %>% pull(value_nat)),
  "canopy" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("canopy", coefs) & !grepl("sq", coefs)) %>% 
             pull(value_nat)),
  "canopy_2" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("canopy_sq", coefs)) %>% pull(value_nat)),
  "slope" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("slope", coefs)) %>% pull(value_nat)),
  "herby" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("herby", coefs)) %>% pull(value_nat)),
  "memory" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("memory", coefs) & !grepl("sq", coefs)) %>% 
             pull(value_nat)),
  "memory_2" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("memory_sq", coefs)) %>% pull(value_nat)),
  "sl" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("step", coefs) & !grepl("log", coefs)) %>% 
             pull(value_nat)),
  "log_sl" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("log_step", coefs)) %>% pull(value_nat)),
  "cos_ta" = hour_harmonics_matrix %*% 
    matrix(coefs_clr %>% filter(grepl("turn_a", coefs)) %>% pull(value_nat))
  )

head(harmonics_nat_df_2p)
##   hour     ndvi    ndvi_2   canopy  canopy_2      slope     herby   memory        memory_2           sl     log_sl
## 1  0.0 9.530534 -16.14346 2.099350 -4.041525 -0.2806850 0.1123308 24064517 -41415792435287 -0.002559160 0.04194061
## 2  0.1 9.570825 -16.20680 2.088156 -4.042115 -0.2843202 0.1121408 23612725 -40149347051813 -0.002672645 0.03863821
## 3  0.2 9.599658 -16.25926 2.074921 -4.039583 -0.2877461 0.1119922 23153804 -38855261640541 -0.002782553 0.03534126
## 4  0.3 9.617051 -16.30080 2.059711 -4.033965 -0.2909546 0.1118821 22688790 -37536794743147 -0.002888600 0.03205669
## 5  0.4 9.623059 -16.33146 2.042595 -4.025302 -0.2939379 0.1118074 22218746 -36197304753928 -0.002990512 0.02879152
## 6  0.5 9.617775 -16.35127 2.023650 -4.013644 -0.2966890 0.1117645 21744755 -34840240889741 -0.003088026 0.02555280
##       cos_ta
## 1 -0.2087064
## 2 -0.2300259
## 3 -0.2510494
## 4 -0.2717232
## 5 -0.2919940
## 6 -0.3118097

# turning into a long data frame for plotting with ggplot
harmonics_nat_df_2p_long <- pivot_longer(harmonics_nat_df_2p, cols = !1, names_to = "coef")

Plot the natural scale coefficients. These are now on quite different scales as the covariates are on different scales - particularly the previous space use density, which we’ll omit from this plot.


ggplot() +
  geom_path(data = harmonics_nat_df_2p_long %>% 
              filter(!coef %in% c("memory", "memory_2")),
            aes(x = hour, y = value, colour = coef)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_y_continuous(expression(Time-varying~parameter~values~beta)) +
  scale_x_continuous("Hour", breaks = seq(0,24,2)) +
  scale_color_discrete("Estimate") +
  theme_classic() +
  theme(legend.position = "bottom")

To update the Gamma and von Mises distribution from the tentative distributions from Fieberg et al. (2021), Appendix C, we just do the calculation at each time point.


tentative_shape <- 0.438167
tentative_scale <- 534.3507
tentative_kappa <- 0.1848126

hour_coefs_nat_df_2p <- harmonics_nat_df_2p %>% mutate(shape = tentative_shape + log_sl,
                                                 scale = 1/((1/tentative_scale) - sl),
                                                 kappa = tentative_kappa + cos_ta)

# turning into a long data frame
hour_coefs_nat_long_2p <- pivot_longer(hour_coefs_nat_df_2p, cols = !1, names_to = "coef")

Plot the updated movement parameters - notice that we’re scaling the scale parameter by 1000 to make it more interpretable.


ggplot() +
    geom_path(data = hour_coefs_nat_long_2p %>% 
              filter(coef %in% c("shape", "kappa")),
              aes(x = hour, y = value, colour = coef)) +
  geom_path(data = hour_coefs_nat_long_2p %>%
              filter(coef == "scale"),
              aes(x = hour, y = value/1000, colour = coef)) +
    geom_hline(yintercept = 0, linetype = "dashed") +
    scale_y_continuous("Value of parameter") +
  scale_x_continuous("Hour", breaks = seq(0,24,2)) +
    scale_color_discrete("Estimate") +
    theme_classic() +
    theme(legend.position = "right")

Checking temporally dynamic movement parameters against observed data

Here we sample from the updated step length distribution (we could also follow the same process for turning angles) to generate a distribution for each hour of the day, and assess how well it matches the observed step lengths.


# summarise the observed step lengths
movement_summary <- buffalo_data %>% filter(y == 1) %>%  group_by(id, hour) %>% 
  summarise(mean_sl = mean(sl), median_sl = median(sl))
## `summarise()` has grouped output by 'id'. You can override using the `.groups` argument.

# number of samples from the Gamma distribution
n <- 1e5

# create some empty variables to store the results
gamma_dist_list <- vector(mode = "list", length = nrow(hour_coefs_nat_df_2p))
gamma_mean <- c()
gamma_median <- c()
gamma_ratio <- c()

for(hour_no in 1:nrow(hour_coefs_nat_df_2p)) {
  
  # sample from the Gamma distribution
  gamma_dist_list[[hour_no]] <- rgamma(n, shape = hour_coefs_nat_df_2p$shape[hour_no], 
                                       scale = hour_coefs_nat_df_2p$scale[hour_no])
  
  # summarise
  gamma_mean[hour_no] <- mean(gamma_dist_list[[hour_no]])
  gamma_median[hour_no] <- median(gamma_dist_list[[hour_no]])
  gamma_ratio[hour_no] <- gamma_mean[hour_no] / gamma_median[hour_no]

}

gamma_df_2p <- data.frame(model = "2p", 
                          hour = hour_coefs_nat_df_2p$hour, 
                          mean = gamma_mean, 
                          median = gamma_median, 
                          ratio = gamma_ratio) 

mean_sl_2p <- ggplot() +
  geom_path(data = movement_summary, 
            aes(x = hour, y = mean_sl, colour = factor(id))) +
  geom_path(data = gamma_df_2p, aes(x = hour, y = mean), 
            colour = "red", linetype = "dashed") +
  scale_x_continuous("Hour", breaks = seq(0,24,2)) +
  scale_y_continuous("Mean step length") +
  scale_colour_viridis_d("Buffalo") +
  ggtitle("Observed and modelled mean step length",
          subtitle = "Two pairs of harmonics") +
  theme_classic() +
  theme(legend.position = "none")

mean_sl_2p


median_sl_2p <- ggplot() +
  geom_path(data = movement_summary, 
            aes(x = hour, y = median_sl, colour = factor(id))) +
  geom_path(data = gamma_df_2p, aes(x = hour, y = median), 
            colour = "red", linetype = "dashed") +
  scale_x_continuous("Hour", breaks = seq(0,24,2)) +
  scale_y_continuous("Median step length") +
  scale_colour_viridis_d("Buffalo") +
  ggtitle("Observed and modelled median step length",
          subtitle = "Two pairs of harmonics") +
  theme_classic() +
  theme(legend.position = "none")

median_sl_2p

Creating selection surfaces

As we have both quadratic and harmonic terms in the model, we can reconstruct a ‘selection surface’ to visualise how the coefficient changes through time.

To illustrate, we have a coefficient for the linear term and a coefficient for the quadratic term for every hour of the day (or for every time point that we used to reconstruct the temporal dynamics, which is 0.1 hour in this case). Using these coefficients, we can plot the the quadratic selection curve at the scale of the environmental variable.

Using the natural scale coefficients of NDVI from the model, let’s take the coefficient from Hour 3 and plot the quadratic curve


# first get a sequence of NDVI values, starting from the minimum observed in the data to the maximum
ndvi_min <- min(buffalo_data$ndvi_temporal, na.rm = TRUE)
ndvi_max <- max(buffalo_data$ndvi_temporal, na.rm = TRUE)
ndvi_seq <- seq(ndvi_min, ndvi_max, by = 0.01)
# ndvi_seq

Similar to the harmonics, we add the linear and quadratic terms to get the quadratic curve


coef_hour <- 3

# we can separate to the linear term
ndvi_linear_coef <- hour_coefs_nat_df_2p$ndvi[which(hour_coefs_nat_df_2p$hour == coef_hour)]
ndvi_linear_coef
## [1] 6.748646

# multiply the coefficient by the values of NDVI 
# the coefficient is positive so this will be a positive relationship
ndvi_linear_selection <- ndvi_linear_coef * ndvi_seq 
plot(x = ndvi_seq, y = ndvi_linear_selection,
     main = "Selection for NDVI - linear term at hour 3",
     xlab = "NDVI", ylab = "Estimated selection")
lines(ndvi_seq, rep(0,length(ndvi_seq)), lty = "dashed")



# and the quadratic term - this term is negative so it will create a 'hat' rather 
# than a 'bowl', with decreasing selection at the extremes
ndvi_quadratic_coef <- hour_coefs_nat_df_2p$ndvi_2[which(hour_coefs_nat_df_2p$hour == coef_hour)]
ndvi_quadratic_coef
## [1] -14.03231

# multiply the quadratic coefficient by the values of NDVI^2
ndvi_quadratic_selection <- ndvi_quadratic_coef * (ndvi_seq ^ 2)
plot(x = ndvi_seq, y = ndvi_quadratic_selection,
     main = "Selection for NDVI - quadratic term at hour 3",
     xlab = "NDVI", ylab = "Estimated selection")
lines(ndvi_seq, rep(0,length(ndvi_seq)), lty = "dashed")


# and sum them both together
ndvi_sum_selection <- ndvi_linear_selection + ndvi_quadratic_selection
plot(x = ndvi_seq, y = ndvi_sum_selection,
     main = "Selection for NDVI - linear and quadratic terms at hour 3",
     xlab = "NDVI", ylab = "Estimated selection")
lines(ndvi_seq, rep(0,length(ndvi_seq)), lty = "dashed")

We can see that the quadratic curve at Hour 3 shows highest selection for NDVI values slightly above 0.2, and the coefficient is greater than 0 for NDVI values from around 0 to 0.45.

For brevity we won’t plot the linear and quadratic terms separetely, but we can do so if needed.

Now for Hour 12


hour_no <- 12

# we can separate to the linear term
ndvi_linear_selection <- hour_coefs_nat_df_2p$ndvi[which(hour_coefs_nat_df_2p$hour 
                                                         == hour_no)] * ndvi_seq
# plot(x = ndvi_seq, y = ndvi_linear_selection,
#      main = "Selection for NDVI - linear term",
#      xlab = "NDVI", ylab = "Estimated selection")

# and the quadratic term
ndvi_quadratic_selection <- (hour_coefs_nat_df_2p$ndvi_2[which(hour_coefs_nat_df_2p$hour 
                                                               == hour_no)] * (ndvi_seq ^ 2))
# plot(x = ndvi_seq, y = ndvi_quadratic_selection,
#      main = "Selection for NDVI - quadratic term",
#      xlab = "NDVI", ylab = "Estimated selection")

# and the sum of both
ndvi_sum_selection <- ndvi_linear_selection + ndvi_quadratic_selection
plot(x = ndvi_seq, y = ndvi_sum_selection,
     main = "Selection for NDVI - linear and quadratic terms at hour 12",
     xlab = "NDVI", ylab = "Estimated selection")
lines(ndvi_seq, rep(0,length(ndvi_seq)), lty = "dashed")

Whereas for Hour 12, the coefficient shows highest selection for NDVI values around 0.5, and the coefficient is positive for all NDVI values above 0. The magnitude of selection is also greater than for Hour 3.

We can imagine viewing these plots for every hour of the day, where each hour has a different quadratic curve, but this would be a lot of plots. Another way to look at it is as a 3D surface, where the x-axis is the hour of the day, the y-axis is the NDVI value, and the z-axis (colour) is the coefficient value.

All we have to do is index over every time point and use the linear and quadratic coefficient values to reconstruct the quadratic curves, and then combine those curves into a 3D surface.

We use the same process as above, but now we index across every time point.

In the plotting we’ll add some contours for the 0 line, to show where the coefficient is positive and negative (conditional on all other covariates).

Selection surface for NDVI


ndvi_min <- min(buffalo_data$ndvi_temporal, na.rm = TRUE)
ndvi_max <- max(buffalo_data$ndvi_temporal, na.rm = TRUE)
ndvi_seq <- seq(ndvi_min, ndvi_max, length.out = 75)

# Create empty data frame
ndvi_fresponse_df <- data.frame(matrix(ncol = nrow(hour_coefs_nat_df_2p), 
                                       nrow = length(ndvi_seq)))
# index across all time points
for(i in 1:nrow(hour_coefs_nat_df_2p)) {
  # Extract the coefficient values for the linear and quadratic terms and 
  # multiply by the NDVI values
  # Assign the vector as a column to the dataframe
  ndvi_fresponse_df[,i] <- (hour_coefs_nat_df_2p$ndvi[i] * ndvi_seq) + 
    (hour_coefs_nat_df_2p$ndvi_2[i] * (ndvi_seq ^ 2))
}

ndvi_fresponse_df <- data.frame(ndvi_seq, ndvi_fresponse_df)
colnames(ndvi_fresponse_df) <- c("ndvi", hour)
ndvi_fresponse_long <- pivot_longer(ndvi_fresponse_df, cols = !1, names_to = "hour")

ndvi_contour_max <- max(ndvi_fresponse_long$value) # 0.7890195
ndvi_contour_min <- min(ndvi_fresponse_long$value) # -0.7945691
ndvi_contour_increment <- (ndvi_contour_max-ndvi_contour_min)/10

ndvi_quad_2p <- ggplot(data = ndvi_fresponse_long, aes(x = as.numeric(hour), y = ndvi)) +
  geom_point(aes(colour = value)) + # colour = "white"
  geom_contour(aes(z = value), 
               breaks = seq(ndvi_contour_increment, ndvi_contour_max, ndvi_contour_increment), 
               colour = "black", linewidth = 0.25, linetype = "dashed") +
  geom_contour(aes(z = value), 
               breaks = seq(-ndvi_contour_increment, ndvi_contour_min, -ndvi_contour_increment), 
               colour = "red", linewidth = 0.25, linetype = "dashed") +
  geom_contour(aes(z = value), breaks = 0, colour = "black", linewidth = 0.5) +
  scale_x_continuous("Hour", breaks = seq(0,24,6)) +
  scale_y_continuous("NDVI value", breaks = seq(-1, 1, 0.25)) +
  scale_colour_viridis_c("Selection") +
  # ggtitle("Normalised Difference Vegetation Index (NDVI)") +
  theme_classic() +
  theme(legend.position = "right")

ndvi_quad_2p

We can now see how the selection for NDVI changes throughout the day. We can clearly see the diurnal pattern of selection for NDVI, where the buffalo are seeking higher values of NDVI (centred on around 0.5) during the middle of the day, indicated by the positive coefficients, and relatively low selection for lower values of NDVI (centred on around 0.2) during the early morning and late afternoon, and negative coefficients for NDVI values above around 0.4, and strong selection against high values of NDVI. The selection against high NDVI values correlates with the movement dynamics, as they are likely moving through more open areas.

Canopy cover


canopy_min <- min(buffalo_data$canopy_01, na.rm = TRUE)
canopy_max <- max(buffalo_data$canopy_01, na.rm = TRUE)
canopy_seq <- seq(canopy_min, canopy_max, length.out = 75)

# Create empty data frame
canopy_fresponse_df <- data.frame(matrix(ncol = nrow(hour_coefs_nat_df_2p), 
                                         nrow = length(canopy_seq)))
for(i in 1:nrow(hour_coefs_nat_df_2p)) {
  # Assign the vector as a column to the dataframe
  canopy_fresponse_df[,i] <- (hour_coefs_nat_df_2p$canopy[i] * canopy_seq) + 
    (hour_coefs_nat_df_2p$canopy_2[i] * (canopy_seq ^ 2))
}

canopy_fresponse_df <- data.frame(canopy_seq, canopy_fresponse_df)
colnames(canopy_fresponse_df) <- c("canopy", hour)
canopy_fresponse_long <- pivot_longer(canopy_fresponse_df, cols = !1, names_to = "hour")

canopy_contour_min <- min(canopy_fresponse_long$value) # 0
canopy_contour_max <- max(canopy_fresponse_long$value) # 2.181749
canopy_contour_increment <- (canopy_contour_max-canopy_contour_min)/10

canopy_quad_2p <- ggplot(data = canopy_fresponse_long, aes(x = as.numeric(hour), y = canopy)) +
  geom_point(aes(colour = value)) +
  geom_contour(aes(z = value), 
               breaks = seq(canopy_contour_increment, canopy_contour_max, canopy_contour_increment), 
               colour = "black", linewidth = 0.25, linetype = "dashed") +
  geom_contour(aes(z = value),
               breaks = seq(-canopy_contour_increment, canopy_contour_min, -canopy_contour_increment),
               colour = "red", linewidth = 0.25, linetype = "dashed") +
  geom_contour(aes(z = value), breaks = 0, colour = "black", linewidth = 0.5) +
  scale_x_continuous("Hour", breaks = seq(0,24,6)) +
  scale_y_continuous("Canopy cover", breaks = seq(0, 1, 0.25)) +
  scale_colour_viridis_c("Selection") +
  # ggtitle("Canopy Cover") +
  theme_classic() +
  theme(legend.position = "right")

canopy_quad_2p

The selection surface for Canopy cover tells a similar story to NDVI, with the buffalo selecting for higher values of canopy cover during the middle of the day, and lower values during the early morning and late afternoon, suggesting that they are seeking more cover during the hotter parts of the day. Overall the values of the coefficient are not as high as for NDVI, and only range between around -1 and +0.5, whereas for NDVI the coefficient ranges from -4 to +4, suggesting that NDVI is more influential.

Previous space use density


memory_min <- min(buffalo_data$kde_memory_density, na.rm = TRUE)
# memory_max <- max(buffalo_data$kde_memory_density, na.rm = TRUE)

# take the 95% quantile instead of the max as the selection coefficients get quite large
memory_max <- quantile(buffalo_data |> filter(y == 1) |> 
                         pull(kde_memory_density), probs = 0.95, na.rm = TRUE)
memory_seq <- seq(memory_min, memory_max, length.out = 75)

# Create empty data frame
memory_fresponse_df <- data.frame(matrix(ncol = nrow(hour_coefs_nat_df_2p), 
                                         nrow = length(memory_seq)))
for(i in 1:nrow(hour_coefs_nat_df_2p)) {
  # Assign the vector as a column to the dataframe
  memory_fresponse_df[,i] <- (hour_coefs_nat_df_2p$memory[i] * memory_seq) + 
    (hour_coefs_nat_df_2p$memory_2[i] * (memory_seq ^ 2))
}

memory_fresponse_df <- data.frame(memory_seq, memory_fresponse_df)
colnames(memory_fresponse_df) <- c("memory", hour)
memory_fresponse_long <- pivot_longer(memory_fresponse_df, cols = !1, names_to = "hour")

memory_contour_min <- min(memory_fresponse_long$value) # 0
memory_contour_max <- max(memory_fresponse_long$value) # 2.181749
memory_contour_increment <- (memory_contour_max-memory_contour_min)/10

memory_quad_2p <- ggplot(data = memory_fresponse_long, aes(x = as.numeric(hour), y = memory)) +
  geom_point(aes(colour = value)) +
  geom_contour(aes(z = value), 
               breaks = seq(memory_contour_increment, memory_contour_max, 
                            memory_contour_increment), 
               colour = "black", linewidth = 0.25, linetype = "dashed") +
  geom_contour(aes(z = value),
               breaks = seq(-memory_contour_increment, 
                            min(-memory_contour_increment, memory_contour_min), 
                            -memory_contour_increment),
               colour = "red", linewidth = 0.25, linetype = "dashed") +
  geom_contour(aes(z = value), breaks = 0, colour = "black", linewidth = 0.5) +
  scale_x_continuous("Hour", breaks = seq(0,24,6)) +
  scale_y_continuous("Previous space use", labels = scientific) +
  scale_colour_viridis_c("Selection") +
  # ggtitle("Relationship to previous space use") +
  theme_classic() +
  theme(legend.position = "right")

memory_quad_2p
## Warning: `stat_contour()`: Zero contours were generated
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## Warning: `stat_contour()`: Zero contours were generated
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf

The relationship to previous space use is a bit more complicated to interpret, and for other more or less harmonics the trend changes, which is unlike NDVI and canopy cover, which are stable from 1 or more pairs of harmonics. There is consistent selection for higher previous space use density in the evening (between around 6pm and 10pm), suggesting that they are attracted to areas they have previously used, which may be foraging areas or similar, although this would be interesting to investigate in the field. It should also be noted that there is no 0 contour, as the entire selection surface is above 0, indicating that buffalo will, on average, always select for areas that have been previously. This selection for areas of previous space use lend the home ranging behaviour in the simulations. The coefficient values range from near 0 to up around +4, suggesting that previous space use is quite influential, particularly in the evening.

References

Fieberg, John, Johannes Signer, Brian Smith, and Tal Avgar. 2021. “A ’How to’ Guide for Interpreting Parameters in Habitat-Selection Analyses.” The Journal of Animal Ecology 90 (5): 1027–43. https://doi.org/10.1111/1365-2656.13441.

Session info


sessionInfo()
## R version 4.2.1 (2022-06-23 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_New Zealand.utf8  LC_CTYPE=English_New Zealand.utf8    LC_MONETARY=English_New Zealand.utf8
## [4] LC_NUMERIC=C                         LC_TIME=English_New Zealand.utf8    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] cowplot_1.1.1       ggExtra_0.10.1      ggh4x_0.2.6         Rfast_2.0.7         RcppZiggurat_0.1.6 
##  [6] formatR_1.14        scales_1.2.1        glmmTMB_1.1.8       clogitL1_1.5        Rcpp_1.0.10        
## [11] ecospat_3.5         TwoStepCLogit_1.2.5 survival_3.5-5      viridis_0.6.2       viridisLite_0.4.2  
## [16] matrixStats_1.0.0   patchwork_1.1.2     ggpubr_0.6.0        adehabitatHR_0.4.21 adehabitatLT_0.3.27
## [21] CircStats_0.2-6     boot_1.3-28.1       MASS_7.3-59         adehabitatMA_0.3.16 ade4_1.7-22        
## [26] sp_1.6-0            ks_1.14.0           beepr_1.3           tictoc_1.2          terra_1.7-23       
## [31] amt_0.2.1.0         lubridate_1.9.2     forcats_1.0.0       stringr_1.5.0       dplyr_1.1.2        
## [36] purrr_1.0.1         readr_2.1.4         tidyr_1.3.0         tibble_3.2.1        ggplot2_3.4.2      
## [41] tidyverse_2.0.0    
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.3             tidyselect_1.2.0       lme4_1.1-32            htmlwidgets_1.6.2     
##   [5] grid_4.2.1             pROC_1.18.0            munsell_0.5.0          codetools_0.2-19      
##   [9] ragg_1.2.5             units_0.8-1            miniUI_0.1.1.1         withr_2.5.0           
##  [13] audio_0.1-10           colorspace_2.1-0       highr_0.10             knitr_1.42            
##  [17] rstudioapi_0.14        ggsignif_0.6.4         Rdpack_2.4             labeling_0.4.2        
##  [21] emmeans_1.8.5          TeachingDemos_2.12     bit64_4.0.5            farver_2.1.1          
##  [25] coda_0.19-4            vctrs_0.6.2            generics_0.1.3         TH.data_1.1-2         
##  [29] circular_0.4-95        xfun_0.39              timechange_0.2.0       randomForest_4.7-1.1  
##  [33] R6_2.5.1               isoband_0.2.7          cachem_1.0.7           reshape_0.8.9         
##  [37] promises_1.2.0.1       vroom_1.6.1            multcomp_1.4-23        nnet_7.3-18           
##  [41] gtable_0.3.3           mda_0.5-3              sandwich_3.0-2         rlang_1.1.0           
##  [45] systemfonts_1.0.4      splines_4.2.1          rstatix_0.7.2          TMB_1.9.10            
##  [49] earth_5.3.2            broom_1.0.4            checkmate_2.1.0        biomod2_4.2-2         
##  [53] yaml_2.3.7             reshape2_1.4.4         abind_1.4-5            backports_1.4.1       
##  [57] httpuv_1.6.9           Hmisc_5.0-1            tools_4.2.1            nabor_0.5.0           
##  [61] ellipsis_0.3.2         raster_3.6-20          jquerylib_0.1.4        RColorBrewer_1.1-3    
##  [65] proxy_0.4-27           plyr_1.8.8             base64enc_0.1-3        classInt_0.4-9        
##  [69] rpart_4.1.19           zoo_1.8-12             cluster_2.1.4          tinytex_0.48          
##  [73] magrittr_2.0.3         data.table_1.14.8      mvtnorm_1.1-3          fitdistrplus_1.1-8    
##  [77] mime_0.12              hms_1.1.3              evaluate_0.20          xtable_1.8-4          
##  [81] mclust_6.0.0           gridExtra_2.3          compiler_4.2.1         KernSmooth_2.23-20    
##  [85] crayon_1.5.2           minqa_1.2.5            htmltools_0.5.5        later_1.3.0           
##  [89] mgcv_1.8-42            tzdb_0.3.0             Formula_1.2-5          DBI_1.1.3             
##  [93] sf_1.0-12              Matrix_1.6-5           car_3.1-2              permute_0.9-7         
##  [97] cli_3.6.1              rbibutils_2.2.13       parallel_4.2.1         pkgconfig_2.0.3       
## [101] numDeriv_2016.8-1.1    foreign_0.8-84         foreach_1.5.2          bslib_0.4.2           
## [105] estimability_1.4.1     plotmo_3.6.2           digest_0.6.31          pracma_2.4.2          
## [109] vegan_2.6-4            rmarkdown_2.21         htmlTable_2.4.1        PresenceAbsence_1.1.11
## [113] shiny_1.7.4            gtools_3.9.4           nloptr_2.0.3           lifecycle_1.0.3       
## [117] nlme_3.1-162           jsonlite_1.8.4         carData_3.0-5          fansi_1.0.4           
## [121] pillar_1.9.0           lattice_0.21-8         fastmap_1.1.1          plotrix_3.8-2         
## [125] glue_1.6.2             gbm_2.1.8.1            iterators_1.0.14       bit_4.0.5             
## [129] class_7.3-21           stringi_1.7.12         sass_0.4.5             maxnet_0.1.4          
## [133] textshaping_0.3.6      poibin_1.5             e1071_1.7-13           ape_5.7-1